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^^ I Abstract. The majority of existing probabilistic model checking case 

^^ . studies are based on well understood theoretical models and distribu- 

^-^' tions. However, real-life probabilistic systems usually involve distribution 

C^ , parameters whose values are obtained by empirical measurements and 

thus are subject to small perturbations. In this paper, we consider pertur- 
bation analysis of reachability in the parametric models of these systems 
m , (i-e-, parametric Markov chains) equipped with the norm of absolute 

distance. Our main contribution is a method to compute the asympto- 
tic bounds in the form of condition numbers for constrained reachability 
probabilities against perturbations of the distribution parameters of the 
system. The adequacy of the method is demonstrated through experi- 
ments. 



1 Introduction 



(N 
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■^^ , Probabilistic model checking is a verification technique that has matured over 

the past decade, and one of the most widely known and used probabilistic model 
^O , checking tools is PRISM [T] . The majority of the reported case studies of prob- 

C ; ' abilistic model checking, including those performed in PRISM, involve systems 

■"sj" , whose stochastic nature arises from well understood theoretical probabilistic dis- 

^D ' tributions, such as the use of a fair coin toss to introduce randomization into an 

algorithm, or the uniform distribution of randomly chosen IP addresses in the 
Zeroconf protocol. More complex, realistic systems, on the other hand, involve 
behaviors or other system characteristics generated by empirical distributions 
that must be encoded via empirically observed parameters. In many cases, these 
^ , distribution parameters are based on finite numbers of samples and are statis- 

Cy_' tical estimations that are subject to further adjustments. Also, the stochastic 

nature of the model (e.g., the failure rate of some hardware component) may be 
varying over time (e.g., the age of the component). The conventional techniques 
and tools of probabilistic model checking, including PRISM, do not provide suf- 
ficient account for systems with distribution parameters. Consider, for instance, 
the setting of automata-based model checking: Given a (probabilistic) model A4 
and an LTL formula ip, the model checker returns a satisfaction probability p 
oi if in M.. However, A4 is just an idealized model of the probabilistic system 
under consideration, and because the real distribution(s) of its parameter(s) may 
be slightly different from those specified in ^Vf, p is merely a referential result 
whilst the actual result may deviate from p to a small but non-trivial extent. 
We elaborate this pitfall in the following two concrete examples. 




Fig. 1. Zeroconf protocol with noisy channels 



Motivating examples. We first present an IPv4 Zeroconf protocol model for a 
network with noisy communication channels. Figure[T]present the protocol model 
that uses a maximum of four 'ok' message probes. Let a be the probability that 
the new host chooses an IP address that has been used already, and x be the 
probability that the probe or its reply is lost due to channel noise or some other 
reason (if any). If an IP address is randomly chosen, then a is equal to m/n, 
where n = 60, 534 is the size of IP address space as specified in Zeroconf and m 
is number of the connected hosts. By contrast, x relies on an ad hoc statistical 
estimation of the lose rate of messages tested in experiments. In reality, it is 
less meaningful to specify a single, constant value of x, as the measurement 
could be affected by factors such as net load, environment noise, temperature, 
measurement time, ect. Instead, x may be given as the expression a;o±a;/i, where 
Xo is the mean value of the measured results and xa specifies the maximal 
perturbation. It is therefore reasonable to express the probability that an IP 
collision happens as the form oi y = yo iyA, where j/o is a referential value for 
the result and y^ specifies the range of perturbation of y. However, although 
the standard model checking techniques allow one to obtain yo by inputting xq, 
they provide little account for the relationship between y^ and xa- 

Another example is an variant of the classic hopping frog example. We sup- 
pose a frog hopping between four rocks, and the (z,j)-entry in the following 
parametric transition matrix (a matrix whose entries are either numbers or vari- 
ables) 
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present the concrete or abstract probability of frog's movement from the ith 
rock to the jth rock. In this matrix, the tuple of variables {zi, Z2, z^, Z4) forms 
a distribution parameter (i.e., satisfying the constraint that Zi > for each 
1 < i < 4 and 21 -f Z2 + 23 + Z4 = 1) and is restricted by the following condition: 
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Table 1. Reachability checking in MCs and PMCs 



Model Input Output 

A4 5*7 U S\ p (idealized reachability probability) 

p (referential reachability probability) 
A4, b-rUbt and Ki (condition numbers) 



where Z\ is a sufficiently small positive number. Equation ([T|) implies that (i) 
the "idealized" value of (zi, 2:2, 23, Z4), which we call a reference to the distribu- 
tion parameter, is i^j'^jjjj) and (ii) a small perturbation of {zi,Z2,Z3, Z4) is 
allowed and measured. A typical model checking problem for this example can 
be stated as "what is the probability that the frog eventually reaches the fourth 
rock without landing on the third rock?" Again, well-established model checking 
techniques do not directly allow one to obtain a solution to this question. 



Approach. In a nutshell, the aforementioned two examples demonstrate that a 
satisfactory model checking result for a probabilistic system with one or more 
perturbed distribution parameters should shed light on the sensitivity of the 
result to the distribution parameters. In this paper, we provide a method to 
compute the (local) asymptotic bounds of the results in terms of condition num- 
bers for reachability checking of probabilistic systems under perturbations. We 
model the probabilistic systems in discrete-time Markov chains (MCsj^ with dis- 
tribution parameters, which are coined as parametric Markov chains (PMCs), 
and introduce the norm of absolute distance to measure the deviation distances 
of their distribution parameters (as exemplified by equation ([T])). The reachabil- 
ity checking is formalized as follows: given a PMC A^* with state space S and 
two sets of states S-?, 5*1 C S*, a reachability problem is phrased as the probabil- 
ity of "reaching states in S\ only via states in S"?" . By adopting notations from 
temporal logic, the problem is denoted by S-fU S\, where U refers to the "until" 
operatoro Two instances of the reachability problem class are mentioned in the 
two motivating examples above. The output of the reachability checking con- 
sists of a referential probabilistic result p G [0, 1] and a condition number k^ G K 
where i G /, an index set, for each distribution parameter. The significance of 
the output is that, if a sufficiently small perturbation Ai, measured by the norm 
of absolute distance, occurs on the parameter whose condition number is Ki for 
each i (^ I, then the actual result is asymptotically bounded by p ± J2iei ^i^i- 
A brief comparison of the reachability checking in MCs and PMCs in terms of 
input and output is presented in Table [TJ 



^ Throughout the paper, unless mentioned otherwise, by MCs we mean discrete-time 
Markov chains. 

^ In fact, the formulation of reachability in the present paper is sightly more gen- 
eral than the standard definition of reachability and sometimes is called constrained 
reachability, since the S-? in S^ U S\ plays a constraining role. 



Our work is related to a line of research |2|3|4|5j . which pursues the per- 
turbation bounds of stationary distributions for MCs. However, to the best 
of our knowledge, this paper is the first paper devoted to the application of 
concepts and methods from perturbation theory to the analysis of quantitative 
verification. To further explain our method, it is beneficial to compare it with 
a standard method for measurement error estimation based on differentiation 
and linear approximation. Suppose a sphere (such as a prototype of balls pro- 
duced by a sporting goods factory) is measured and its radius is 21cm with 
a possible small error within O.OScm. The dependence of the sphere volume 
on the radius is given by y = ^irr^. The problem is to compute volume er- 
ror Va given the radius error ta- We recall a classic method for this problem: 
First, the differential function of F on r is given by dV = inr'^dr. Second, let 
dr = rA = 0.05c?7i (which is significantly small compared with r = 21cm) and 
we obtain Va ^ dV = Att x 21^ x 0.05 « 277cm^. The sensitivity of Va to rA 
is approximated by the ratio ^ = Airr^ « 5, 542 and this expression is useful if 
the value of va is unknown in advance. We aim to develop a similar methodology 
to estimate the perturbations of reachability in PMCs, which is comparable to 
the use of differentiation and linear approximation in estimating the error of the 
ball volume described above. 



Organization. The remainder of the paper is organized as follows. In the next 
section (Section[2]), we present the formulations of PMCs and introduce the norm 
of absolute distance for probabilistic distributions. For presentation purposes, we 
separate the treatment of PMCs into that of basic PMCs, which have a single 
distribution parameter, and non-basic PMCs, which have multiple distribution 
parameters. In Section [3l we deal with basic PMCs by establishing a method 
for computing their asymptotic bounds, in particular, condition numbers for the 
given reachability problems. In SectionlH we generalize the computation method 
to handle non-basic PMCs . We evaluate our approach by experiments in Section 
[5] Related work is discussed in Section [6] In Section [3 the paper is concluded 
and several future research directions are outlined. 



2 Parametric Markov Chains 

In this section, we define the formal models of PMCs, which are parametric vari- 
ants of MCs. Informally speaking, a PMC is obtained from an MC by replacing 
the non-zero entries in one or more rows of its probabilistic transition matrix by 
mutually distinct variables. 

Let X be a symbolic vector of variables, vector variable for short, (xi , . . . ,Xk) 
for some k > 2 such that Xj-^ ^ Xj^ if ji ^ J2- Similar to numerical vectors, we 
use x[j] to denote the variable in the ith entry of x. Two vector variables are 
disjoint if they share no common variables. Let (xi)ig/ be a sequence of pair-wise 
disjoint vector variables for some index set 1^0. Let x* be an extension of x 
obtained by inserting the number zero into x, i.e.. 



x* = (0,...,0,x[l],...,x[j],0,...,0,x[j + l],...,x[fc],0,...,0) , 

ioO's ijO's /fcO's 

for some Iq,- ■ ■ ,lk G N. Such x* is called a distribution parameter. We extend 
each vector variable in (xi)ig/ in similar ways to obtain a sequence of distribution 
parameters such that |x* | = |x* | > max(/) for each ii,i2 G I (note that 
|xij and |xi2 1 are not necessarily the same). Let c S (0,1)'^ be a probability 
distribution, where k = |x|. Thus, X^i^i ^[j] = 1- Let c* be an 'extension' of c 
in the same way as x* of x, i.e., 

c* = (0,...,0,c[l],...,c[i],0,...,0,c[j + l],...,c[fc],0,...,0) . 

ioO's IjO's IfcO's 

We call c* a reference to x*. Similarly, let each vector in (c*)ig/ be a reference to 
the corresponding vector variable in (x*)ig/. For convenience, we write (xi)ig/ 
(resp., (x*)ig/) as x/ (resp., x}), and {ci)iei (resp., {c*)iei) as c/ (resp., c}). 
We also call c} a reference to x/. 

Let P be a transition matrix such that V[i,j] — c*[j] for each i e I. A para- 
m,etric transition matrix 'P(xj) of P is obtained by replacing c* (the probability 
distribution in the ith row of V) with x* for each i G I. We call P a reference 
matrix to V{x*j). If the distribution parameters are clear in the context, we may 
just write V^ for short. 

Definition 1 A parametric Markov chain (PMC) is given by the tuple 

where i is an initial distribution, 'P(x|) is a parametric transition matrix, and 
c*j is a reference to xj . 

Although c/ can be derived from V and x|, for readability we present it 
explicitly in the formulation of 7W*. The underlying MC of Ai^, is 7W = i'-i'P)- 
We do not specify the state space 5 for A^* and Ai. But throughout the paper, 
we assume that S = {1, . . . , \t\}. 

As promised earlier, we introduce a statistical distance measurement between 
distribution parameters and their references, which is given by the norm of ab- 
solute distance (also called total variation). 

Definition 2 The statistical distance for yVf* is given by \\ ■ \\ such that ||a|| = 
X]i=i I^WI /o'^ 'J^l/ vector a. 

By definitions, the scalar function ||x* — c* || is the same as the scalar function 
||x — c||. If / is a singleton, we also call the PMC a basic PMC and denote it 
as A^t- In other words, a basic PMC TW^ is a PMC with a single distribution 
parameter. 



We now present examples of PMCs. The first example of PMCs A^/ is for the 
hopping frog. Its parametric transition matrix (i.e. a 4 x 4 symbolic matrix) of 
which is presented in the Introduction and we denote it as V^ . In 'P*^, the tuple 
(zi, ^2,-23, Z4) is given as the only distribution parameter. We let the reference to 
the parameter be (0.375, 0.125, 0.25, 0.25). In words, ideally, the probabilities for 
the frog to jump from the first rock to the first, second, third, and fourth rocks 
are 0.375, 0.125, 0.25, and 0.25, respectively. Such a PMC is denoted by A^?. 
Additionally, we let the initial distribution in A4? be t^s = (0.25, 0.25, 0.25, 0.25), 
which means that all rock have an equal probability to be the frog's initial 
position. Clearly, J^} is a basic PMC. 

Another example of PMCs M'f is for the noisy Zeroconf. For illustration 
purposes, a probabilistic transition system with a parameter x is provided in 
Figure [TJ The formulation of M^J according to Definition [1] deviates from the 
transition system because of the use of distribution variables. Following the 
definition, we let the sequence of distribution parameters be {xi^x'j)f^i^ where 
Xi + x'^ = 1 for each 1 < i < 4. The parametric transition matrix V^^ is given by 
the following 7x7 symbolic matrix: 
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where the (z, j)-entry in the matrix gives the transition probability from the ith 
state to the jth state. Because the intended initial state is the first state, the 
initial distribution l^^ is (1, 0, . . . , 0). The references are obtained by placing 0.75 
and 0.25 in the positions of Xi and x^, respectively, for each 1 < i < 4. In other 
words, we suppose that under idealized conditions the chances of not receiving 
a reply in four probes are equivalently 0.25. 

3 Perturbation Analysis of Basic PMCs 

From this section, we commence the perturbation analysis of reachability prob- 
lems in PMCs. For presentation purposes, in this section we deal with basic 
PMCs. Recall that a basic PMC has a single distribution parameter. Our main 
goal is to establish a method to compute an asymptotic bound, in particular, a 
condition number for a given reachability problem in a basic PMC against the 
perturbation of its sole distribution parameter. In the next section, we generalize 
the method to the setting of PMCs in general. 



3.1 Perturbation Function 

The reachability problem S-fU S\ in A4^ with state space S such that 5? US'! ^5* 



is formulated in the same way as in A4. For convenience, we let S"? = {1, 



. . ,n7 



?} 
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and S\ = {ni, . . . , \l\}, where < riv < ni < |i|. Thus, S-? n S\ = 0II We cah S? 
the constraint set of S->\S\ and 5: its destination set. 

We use X* to denote the sub-tuple of x* , the only distribution parameter in 
A4^^, such that x* consists of the first n? entries in x*, and x* to denote the 
expression x*[n!] + . . . + x*[|t|]. Similarly, let c* be the the first n? entries of c 
and c* be the sum c*[n!] + . . . + c*[|t|]. Also, let l-? be the first n? entries of t. 

Let p be a vector such that |p| = ti? and p[i] is the satisfaction probability 
of S-fU S{ in state i, for each i G S-/. Informally, let 

717 |t| 

pW-E^[*'J']-p[^1+ E^[*'^"] ' (2) 

for each 1 < z < n?. We rewrite the equation system given in ^ as 

p = A • p + b , (3) 

where A is an nl x n? matrix such that A[i,j] — V[i,j], and b is an n?- vector 

such that h[i] = J2p=n, 'P[hj']- We define p as the least fixed point satisfying 
equation ([3]). 

Lemma 3 p is computed by 

CXD 

p = EA^-b. (4) 

In the following, we define parametric variants of A and b specified in equa- 
tions ^ and Q: Let Ax» be the parametric matrix obtained by replacing the 
first row in A with the tuple x.* , and bj* be the parametric vector by replacing 
the element in the first entry of b with the polynomial x*. (Ax» can be seen as 
the up-left n-/ x n? sub-matrix of 'P(x*).) The selection of the first row (resp. en- 
try) in A (resp. b) should not harm the generality, because we can freely change 
the row positions in V in advance. If x*, n-r and m are clear in the context, we 
also write Ax.* and bj. as A* and b*, respectively. 

As an example, consider the following model checking problem of the hop- 
ping frog (which has already been mentioned in Sections 1 and 2): What is the 
probability of reaching the fourth rock without landing on the third one? In this 
problem, the constraint set is {1,2} and the destination set is {4}. Recall that 
the only distribution parameter in A4^ is (01,^2,2:3,24). Thus, the parametric 
matrix and the parametric vector are respectively given by 

\ 8 8 / V 4 

Let V = [0, l]*^ where fc = |x| and U = {a e V | Y.7=i aW = !}• 



^ It should be noted that this assumption does not impose any theoretical restriction 
on the reachability problem. If S? n Si 7^ 0, we can carry out the analysis based on 
(S^\Si)USi. 



Definition 4 The perturbation function of x (for the problem StU S\ in M.^ 
with reference c) is defined by p :\ ^f [—1,1] such that 

n? oo 

oo 

= .vE(Ai.rbxr-A^"-b) . (5) 

The perturbation function p captures tlie effect of any small variation of x 
with respect to c on the satisfaction probability of the problem S->IA S\ in A4^. 
For convenience, we call c the reference of p. 

3.2 Asymptotic Bounds 

There are various ways to express the asymptotic bounds. We adopt the most 
basic way: The bounds are given by the so-called absolute condition numbers. 
In Section [B] we briefly discuss the terminologies of perturbation bounds and 
condition numbers. 

Let A > represent the perturbation distance of a distribution parameter. 
In reality, we usually assume Z\ as a sufficiently small positive number. The 
following auxiliary definition captures the variation range of p with respect to 
the perturbation distance A of the distribution parameter x* . 

Definition 5 The variation range of p (with reference c) by A is the set 

p(zi) = {p(a) I ||a-c||<AaeU} . (6) 

It is not hard to see that 'p{A) is an interval. The existence of a condition 
number for p depends on the differentiability of p. The following proposition 
confirms that p enjoys this property in a "neighborhood" of c. Recall that we 
have assumed |x| = k. 

Proposition 6 p is differentiable at c, i.e. 

p(y + c) = h . y + e{y) , (7) 

for some h G M'' and 6 : R'' ^ R such that lim||y[|^o ^(y)/l|y|l = 0. 

In other words, h • y is used as the linear approximation of p at a point 
sufficiently close to c, and we write p{y + c) « h-y. Let max(h) — max{h[i] | 1 < 
i < \h\} and min(h) = min{h[z] | 1 < i < |h|}. 

Theorem 7 The asymptotic bound of p is given by the condition number 

K = lim supjl I a; e p((5),0 < 5 < zij . (8) 

Then, the number k exists and, moreover, 

K = -(max(h) — min(h)) , (9) 

where h is given in Equation ([7]). 



According to the definition of k in Theorem [7] (in particular, equation ([S])). 
mathematically, if the parameter x in a basic PMC A^^ with reference c varies 
an infinitesimally small A from c in terms of the absolute distance, then the 
perturbation of the reachability checking result, p{A), is estimated to be within 
±kZ\, where k is the condition number of p. We test the applicability of such k 
in experiments in Section 5. 

The definition of k captures the sensitivity of p to x: How does p change if we 
perturb x? A closely related problem is phrased as this: How much do we have to 
perturb x to obtain an approximation of g — in other words, what is the backward 
error of p? The following proposition gives a "backward" characterization of the 
asymptotic bound n, which, by its formulation, pursues the infimum of variations 
of X (or equivalently, the supremum of their reciprocals) that can cause the given 
perturbation of p. 

Proposition 8 Let k be defined as in Theorem^ 

K = lim sup \^\0<y<x,ye p{S) \ . (10) 

a;-^0 L ' J 

In the following, we present a method to compute the linear approximation 
of p. We write X^i^o ^^ ^® S ^- We define two specific parametric components: 
Dy* is the n? x n? parametric matrix such that Dy* [1, j] = y*[j] and all other 
entries in it are zero; and dp* is the n? parametric vector such that dp* [1] = y* 
and all other entries in it are zero. 

Theorem 9 Let 

Q = 5]A.Dy..^A.b + ^A.dy. . (11) 

Then, p(y + c) w t? • Q. 

Theorems [7] and |9] together provide algorithmic techniques for computing the 
condition number n for A^^ and the reachability problem S->IAS\. 

4 Perturbation Analysis of General PMCs 

In this section, we generalize the method developed in the previous section from 
basic PMCs to general PMCs that may have multiple distribution parameters. 
For general PMCs, perturbations of the parameters may vary either proportion- 
ally or independently, yielding two forms of asymptotic bounds or, equivalently, 
condition numbers. However, it turns out that the two forms of bounds coincide. 

4.1 Directional Conditioning 

For general PMCs, we need to handle multiple distribution parameters. The 
reachability problem SiU S\ in A4* is the same as in 7W^. In this subsection, we 
suppose their perturbations are subject to a prescribed ratio, i.e., proportionally. 
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Hence, we associate a function w : I ^>- [0,1] of to x/ such that J2iei '^(*) ~ -'^■ 
Such w is called a direction of x/ . 

To enable the formal treatment, we first give some notation definitions. For 
each i e I, let V, = [0, l]*^' where h = |xi|. Let U, = {a e V, | X^jLi »[«] = !}■ 
li I = {ii, . . . , im}, then V/ denotes the cartesian space V^^ x . . . x V^^^ and U/ 
has the similar meaning. Let n9,n\ have the same meaning as in the previous 
section. We use x* , to denote the sub-tuple of x* which consists of the first n? 
variables in x*, and x*, to denote the expression x*[n!] + ... + x*[|i|]. Let c*, 
be the the first n? entries of C; and c*, be the sum c*[n\] + . . . + c*[|i|]. Let 



,h 



€1, X} 



(x*,) 



(c*,),g/, and cf, = (c*,)ie/. Let A^ 



be the parametric matrix obtained by replacing the ith row with the tuple x* , 
for each i G /, and bx* , be the parametric vector obtained by replacing the 
«th entry with the expression x* , for each i G /. As before, if x* is clear in the 
context, we write A.^ and b.^ for Ax* ,, and bj. , respectively. 

We illustrate these definitions by the example of noisy Zeroconf, whose model 
is a non-basic PMC. Clearly, the pursuit of the problem "what is probability of 
an IP collision?" is equivalent to the problem "what is probability to avoid an 
IP collision?" In the second problem, the constraint set is {1,...,5} and the 
destination set is {7}. The sequence of parameters is {xi, 1 — Xi)^^^^. Thus, 
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The following definition generalizes Definition |4l 

Definition 10 The perturbation function o/x/ (for the problem S-rU S\ in M.^ 
with reference cj) is defined by g : V/ — >■ [—1,1] such that 



g{^j) = £ i[^] g(A^,. ,^ . b;j._, - A^- • b)[*] 

i=l j=0 

00 



(12) 



3=0 



The perturbation function g captures the effect of the small variation of x, 
with respect to c^ for each i e / on the reachability problem S'/U S\ in A^,. We 
also call c/ the reference of g. The following definition generalizes Definition [5] 

Definition 11 The w-direction variation range of g (with stationary vector c) 
against A is the set 



g^{A) = {g{ai) \ ||a, - c,;|| < w{i)A, a, G U„ i G /} 
where a/ = (ai)ig/. 



(13) 
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It is not hard to see that g^,(^) i'' ^^ mterval. Suppose |yi| — \ci\ for each 
i e I and y/ + c/ is the sequence (y^ + c.j)iG/. Let ||y/|| be J2iei Hy^H- Similar 
to p, the following proposition holds for g. 

Proposition 12 g is differentiable at C/, i.e. 

g{yi + cj)^Y.^,-y, + 0'{yj) (14) 

iei 

where h, G R'^ (i e I) and 9' : R^^l^l -^ R such that lim||y^||^o ^'(y/)/l|y/|| = 0. 

We write g{yi -\- c/) w J2iei ^i • Yi and call J2i£i ^i ' y* the linear approxi- 
mation of g at cj. The following theorem generalizes Theorem [71 

Theorem 13 The ly-direction asymptotic bound of g is given by the directional 
condition number 

Ktu = hm sup I - I X e e^ ((5) , < (5 < Z\ ^ . (15) 

Then, n^ exists and, moreover, 

Hw = -^w(i)(max(hi) - min(h,,)) , (16) 

where each h^ is given in (J14p . 

If the distribution parameters vary a small A in the direction w, then we 
can estimate the perturbation of p as ±k^A. In the case that w{i) = 1/|/| for 
each i G I, such k^ is called a uniform condition number. Like the asymptotic 
bounds for basic PMCs, a "backward" characterization of k^, also exists. 

Proposition 14 Let Kw be defined as in Theorem \lS\ . 

Kt„ = lim sup i - I < y < a;, y G '§^{5) \ ■ (17) 

We also provide a method to compute the linear approximation of g. We 
define two specific parametric matrices. Dy» ^ is the n? x n? matrix such that 
Dy» ,^ [j, j] = Ay- ^ [i, j] for each « G / and 1 < j < fc and all other entries in Dy* , 
are zero, dy* ^ is the nv-vector such that dy* ^ [i] — by* ^ [i] for each i d I and 
all other entries in dy* ^ are zero. We have the following generalized theorem of 
Theorem m 

Theorem 15 For each i G I, let 

Q/=^A.Dy;., .J^A-b + J^A-dy.^ . (18) 

Then, g{yi +C7) « t? -Q/. 

Theorems [T3] and [12] together provide an algorithmic method for comput- 
ing the directional condition number «;„, for Ad^, and the reachability problem 



12 

4.2 Parameter-wise Conditioning 

The parameter-wise perturbation analysis handles the independent variations of 
distribution parameters. In this case, to facilitate the perturbation estimation, 
we expect to obtain a condition number for each distribution parameter. It turns 
out that the parameter-wise analysis can be reduced to the directional analysis. 
We use cj[i :— a] denote the sequence of vectors obtained by replacing the 
«th vector in c/ by a. 

Definition 16 The variation range of g projected at Xi against A is the set 

g,{A) = {q{cj[i := a]) \ ||a-c,|| < zi, a e UJ . (19) 

The asymptotic bound of g projected at x^ is given by the condition number 

K, = hm sup|^ I xe ft((5),0<5<Z\| . (20) 

Let Wi be the direction such that Wi{i) = i and Wi{j) = for each j € /\{J}. It 
is easy to see that 'g^ (resp. Ki) is just g^. (resp. k^.). It means that parameter- 
wise bounds are special cases of directional bounds. Moreover, the following 
theorem states that any set of parameter- wise condition numbers conforms to a 
specific directional condition number. 

Tlieorem 17 Let A = J^iei ^i ^''^^ ^(*) = A/ A for each i G I. Then, 

^K,Z\, = K„Zi . (21) 

If the direction of perturbation may not be known in advance, it is more use- 
ful to present the set of parameter- wise bounds. Theorem 1171 provides a mathe- 
matical characterization for parameter-wise perturbations in terms of direction 
perturbations. 

5 Experiments 

We evaluate by experiments, informally speaking, how well the condition num- 
bers handle possible perturbations of the reachability problem for the PMC 
under consideration. Recall that the outcome of the reachability checking algo- 
rithm for a PMC consists of two parts, namely, a referential probabilistic result 
and one or more condition numbers (see Table [ij . The probabilistic result is 
computed by a conventional numerical model checking algorithm. In particular, 
the probability result is given by the equation 

p^ 17 -p + ^ni] + ... + L[\i\] , (22) 

where p is given in equation (|3]). For problems considered in this section, only 
a single condition number will be returned. The condition number is calculated 
by the method presented in the previous sections. 
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Table 2. Experimental data of noisy Zeroconf (xlO^"^) 

T,, , , T,,i-i-, T^-. Condition Variation 

Model Xi Probability Distance Number Ranee 



Mf 


750 


999.024 


- 


Mf 


749 


-.016 


2 


Mt 


752 


+.031 


4 


Mt 


747 


-.048 


6 



7.797 



±.016 
±.031 

±.047 



Table 3. Experimental data of hopping frog (xlO ■^) 

Model Distribution Probability Distance Number Ranse 

X? (375, 125, 250, 250) 

JvC^ (374,124,251,251) 

M^^ (374, 124, 250, 252) 

M^i (377, 125, 248, 250) 

JvCi (377, 125, 250, 248) 

M^i (375, 125, 248, 252) 



M^^ (375, 125, 252, 248) 



500.000 


- 


312.500 


- 





4 


- 


±1.250 


±.623 


4 


- 


±1.250 


±.627 


4 


- 


±1.250 


-.627 


4 


- 


±1.250 


±1.250 


4 


- 


±1.250 


-1.250 


4 


- 


±1.250 



Our experiments proceed as follows, (i) We specify a reachability problem 
for the PMC A^* and compute the referential probabilistic result p and one con- 
dition number k for the problem, although multiple condition numbers may be 
required in other contexts, (ii) By deliberately assigning concrete probability dis- 
tributions to the distribution parameter(s) of A^», we construct several potential 
non-parametric models M.j with sufficiently small statistical distances to A^*. 
(iii) We compute an actual probabilistic result -pj for each hAj and calculate the 
actual distance Aj between the reference(s) in Al, and the corresponding distri- 
bution(s) in M.j. (iv) We compare p — pj, the difference between the referential 
result and an actual result, and ±kAj, the perturbation estimation. 

We perform experiments for the examples of noisy Zeroconf and hopping 
frog, i.e., the PMCs Mf and MI^, in Matlab® [7]. In the hopping frog example, 
for the ease of calculations, we let a in Vf, the parametric transition matrix of 
A4l^ , be 0.2. It is reasonable assume each probe message and its reply are effected 
by the same channel noise level, namely, that the four distribution parameters 
of A41 are perturbed in a uniform direction. Several MCs in both experiments 
are generated by assigning different distributions to the distribution parameters 
in their PMCs. Because the infinite matrix series X^i^o -^"' cannot be computed 
directly, we adopt an approximation by taking the sum of the first hundred items 
in each series encountered. There is no significant truncation error involved the 
numerical calculation in our experiments. We test the reachability problems 
{1, . . . , 5} W {7} in the first experiment (which states "what is the probability to 
avoid an IP collision?" ) and {1,2}U {4} in the second one (which states "what is 
the probability for the frog to reach the fourth rock without landing on the third 
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rock?"). The experimental data are summarized in Figures [2] and |3l respectively. 
In Figure [21 the distance of the perturbed models to the PMC increases. We 
observe that the condition number accounts for the result nicely if the perturbed 
distance is smaller than 0.006. When the distance exceeds 0.006, the perturbation 
of the probabilistic result may drop out of the variation range. In Figure [3l 
several perturbed models with the distance 0.004 to the PMC are inspected. 
The data also demonstrate that the condition number bounds the reachability 
perturbation between a perturbed model and the PMC to a satisfactory degree. 
In particular, the difference between the result for Ai^ (resp., ^A;^) and the 
referential result overlaps with the positive (resp., negative) predicted bound. 

In short, we observe from the experiments that condition numbers adequately, 
although not rigorously, predict the bounds of the reachability checking results 
if distribution parameters perturbations are sufficiently small. 

6 Related Work 

Our work is related to the study of perturbations of MCs, which pursues pertur- 
bation bounds for their stationary distributions, and which can be traced back 
to 1960's. Indeed, from the technical point of view, reachability and the sta- 
tionary distribution have an intrinsic connection, and so do their perturbations. 
Schweitzer [5] gave the first perturbation bound, namely, an absolute condition 
number for the stationary distribution of an MC against its fundamental ma- 
trix (which is defined by the transition matrix of the MC), and motivated a 
variety of subsequent work. Cho and Meyer p] provided an excellent overview 
on various bounds of stationary distributions (all of which are in the forms of 
condition numbers) up to the written time, whilst more recent papers in this 
line of research include [4] and J5] , both of which shed light on new definitions 
and techniques for perturbation bounds. In spite of its relatively long history, 
to the best of our knowledge, the present paper is the first paper that raises the 
perturbation problem in quantitative verification. 

In a more general setting, the main goal of perturbation theory for numerical 
linear algebra jB] is to investigate the sensitivity of given problems with respect 
to the perturbed input data in a prescribed domain and provide various forms 
of perturbation bounds for the problems. One important group of such bounds 
is called asymptotic bounds (also called linear local bounds) , which is further di- 
vided into two subgroups, i.e., absolute condition numbers and relative condition 
numbers. The condition numbers for the reachability verification pursued in the 
present paper belong to the first subgroup. Both subgroups of condition num- 
bers have their significance: absolute condition numbers enjoy a more elegant 
mathematical formulation and are easier to be employed to analyze practical 
problems, whilst relative ones are more important to fioating point arithmetic 
implemented in every computer, which is affected by relative rather than abso- 
lute errors. A detailed classification of these bounds is found in [5] (Chapters 1 
and 2). We leave the analysis of our problem based on other kinds of bounds to 
future work. 
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There have been many efforts that attempted to introduce perturbation er- 
rors to realtime systems, in particular, timed automata. For example, Alur et 
al [10] defined a perturbed semantics for timed automata whose clocks might 
change at some very small rates. They showed that if the automaton has a single 
clock then the language accepted by it under the perturbed semantics is accepted 
by a deterministic automaton under the standard semantics. Bouyer et al [llj 
provided another perturbed semantics, which captures not the perturbations of 
the clock rates but those of the clock constraints. They developed model checking 
techniques of Biichi and LTL properties based on a novel semantic satisfaction 
relation, which is defined via their perturbed semantics and captures — as they 
argued — the intuitive notion "whether the considered property holds for the 
same model implemented in a sufficient (but not infinitely) fast hardware" . Still 
another different problem, which studies realtime systems with general dynamic 
nature, was addressed by Chen et al [T^, who proposed an LTL model check- 
ing method for a continuous-time inhomogeneous (or non-stationary) MC, the 
inhomogeneity of which is predefined by an integrable function. 

Sensitivity of computable problems to their input data are also explored 
in other contexts. Misailovic et al [13! introduced a framework of approximate 
correctness for programs, one utility of which is to analyze the robustness of 
programs with disturbed input, namely, whether small perturbations to the input 
lead to small variations of the output only. Heimdahl et al |14| employed model 
checking techniques in software deviation analysis of control systems. 

7 Conclusions 

Motivated by the pervasive phenomena of perturbations in the modeling and ver- 
ification of real-life probabilistic systems, we study the sensitivity of constrained 
reachability probabilities of those systems, which are modeled by parametric 
variants of discrete-time MCs, to perturbations of their distribution parameters. 
Our contribution is a method to compute the asymptotic bounds in terms of 
absolute condition numbers to characterize the sensitivity. We also conduct ex- 
periments to demonstrate the practical adequacy of the computation method. 

This paper is an initial step of investigating the sensitivity and bounds for 
quantitative verifications of systems under perturbations, and several interesting 
directions for further research are outlined below. First, reachability, in spite of 
its fundamental status in model checking, captures a narrow group of practical 
verification problems (particularly in the probabilistic domain) and, therefore, 
it is desirable to extend the present method to accommodate the general model 
checking problems formalized, for instance, in the LTL formulas. Second, we 
adopt the norm of absolute distance to measure the distance between two prob- 
ability distributions, however, there exist other measurements that are particu- 
larly useful for problems in some specific domains. For example, the well-known 
Kullback-Leibler divergence is widely adopted in information theory [TS]. Fi- 
nally, condition numbers are among several other forms of perturbation bounds. 
An in-depth comparison of their pros and cons is left to future work. 
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Appendix: Proof Details 

Proof (Proof of Lemma\3^. See Theorem 10.15 and its subsequent remark at 
pp.762-764in 17,. 

Proof (Proof of Proposition \^. We refer this proof to the constructive (and 
independent) proof of Theorem |9] below. 

Proof (Proof of Theorem^. Let 

«' = -j^'y^&'^i^) - min(h)) = -(h[ii] - h[z2]) 

for some «i,i2, and our goal is to show that n exists and k = k'. Note that 
|h| = fc. 

First we give a claim. Let q G [—1, l]'^ such that q = and ||q|| ~ 1. We 
claim that 

h • q < k' . 

The proof of the claim is straightforward. 

Let £ > 0. According to Proposition |6l choose Z\ > such that if < ||y|| = 
S < A then 



|h-y-p(y + c)| < 



5e 



Let y = (5q, then 



hq 



p(y + c) 



h • y - p(y + c) 



e 
< - 
- 2 



in particular, 



, p{rei^,i,_+c) 



re. 



< 



where e^ j is the vector such that the zth (resp. jth) entry in eij is 1/2 (resp. —1/2) 
and the other entries (if any) are all zero. Supposing r is sufficiently small, 
rei-^^i^ + c G U. Then, 



Pi'ren,i2 +c) 



e {| I a;ep(5),0<(5<Z\} 



Thus, 



k' - £ < sup II I X e p{S), <s < a\ 
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On the other hand, given < 6' < A, choose Z/i such that ||z|| = 6' , z+c e U, 
and sup(p((5')) < p{z+c)+eS' /2. Without loss of generality, we suppose z = S'q'. 
Thus, 



, ^ p(z + c) e 



Thus, k' + e> ^^^^P}- '' for any 5 such that {)< 5 < A. Hence, 



K' + e > supjl I a; G p((S),0 < J < zij . 
Therefore, we have proved that k, exists and n — k! . 
Proof (Proof of Proposition\B^. We suffice to show the following two propositions: 



1. for each e > there exists a; > such that y/S < k + e whenever < y < a; 
and y £ p((5); 

2. for each e > and a; > there exists < y < a; and S > such that y e 'p{S) 
and K < y/5 + e. 

On the other hand, by the definition of k, 

1 . for each e > there exists A > such that y/S < n + e whenever < 6 < A 
and y G p((5); 

2. for each e > and A > there exists < S < A and y > such that 
y G p(<5) and k < y/6 + e. 

It holds that given any Z\ > we can choose x such that if < y < x and 
y G p((5) then < (5 < Z\; conversely, given any a; > we can choose A such that 
if < ^ < Z\ and y G p(<5) then < y < a;. Therefore, it can be verified that the 
two sets of propositions are equivalent. 



Proof (Proof of Theorem 0). p(x) actually defines a system of non-linear (i.e. 
multi-variable polynomial) series. Given a multi-variable polynomial series, the 
order of a term in the series is the summation of the exponents of all variants 
in it. The smallest order of all terms is called least order of the series. By the 
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definition of /(x), we have that 

p(y + c) 

= t,.(/(y + c)-/(c)) 



^(A + Dy,)^(b + dy,)-^A.b 



Vi=0 



i=0 






V 



^A'.dy.+^A^.Dy,.^A'.b 



V V-i 



+ (^ , 



/ 



where </? is either or a of polynomial whose least order is not smaller than 2. 
Thus, 

lim ---- = . 
llylHo ||y|| 

Let h.- y — L-? ■ {ipi + ■02), and we obtain equation ([7]) in Proposition [HI namely 

p(y + c) « i? • (V'l +'(/'2) ■ 

Proof (Proof of Proposition\W\) . We refer this proof to the proof of Theorem [TSl 
below. 

Proof (Proof of Theorem \13\) . The proof is a generalization of the proof of The- 
orem [71 For any direction w, let 

kL = ^5 E w'(i)(max(h,) - min(h,)) = ^^ E ^'(*)('^»b'»^] " ^^'b?]) > 



iG/ 



iG/ 



for some jl,jf for each i £ I. Our goal is to show that k,^ exists and Kw = k^. 
Note that |hi| = fc for each i. 

We first give a claim. Let q^ G [—1, 1]'^ such that, for each i G /, 'Yl,-i=i ^i[j] = 
and ||qi|| = 1. We claim that 

^w{i){h,-Cii)<K'\I\ . 

The claim is not hard to be verified. 

Let e > 0. According to Proposition[T2l choose A > such that if < ||yi|| = 
6w{i) < Aw{i) for any i £ I, then 



^hj -yj - g{yi + ci) 



< 



Se 



i£l 
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Let Yi = 6w{i)c\i for each i G /. Thus 



^M;(i)(hi-q,) 



Q{yi + c/) 



is/ 



Eie/hj -y,: -(?(y/ + c/) 



< 



in particular, let e/ = (e,ij-2)ig/ where 6,1,2 is the vector such that the jj^th 
(resp. j'fth) entry in it is 1/2 (resp. —1/2) and the other entries (if any) are all 
zero (thus, e,i ,2 is a particular instance of g^); then 

q{5w ■ 6/ + C/) 



e 
< - 
- 2 



Since 5 is sufficiently small, Sw • e/ + c/ G U/. Then 
e((5w • e/ + c/) 



Thus, 



g{| |a;Ge^((5),0<J<Z\} 



K^-e<sup{^ I xG g^^{S),0<S<A} 



On the other hand, given < S' < A, choose z/ = {zi)i^j such that ||zi|| = 
6'w{i) for each i G /, zj + cj G U/, and sup(^((5')) < p{zj + cj)+eS' /2. Without 
loss of generality, we suppose z^ = 6'q'^ for each i G /. Thus, 



h.-q, + £> ^^ + - 



iGl 



Thus, .:,+e>^^%M. Hence, 

K^ +£ > supj^ I a; G e„,((5), < (5 < Z\| . 

Therefore, we have proved that k^ exists and k„ = k^. 

Proof (Proof of Proposition \14^ . The proof is an immediate generalization of the 
proof of Proposition IHl 

Proof (Proof of Theorem \15\) . The proof of this theorem is a generalization of 
that of Theorem ini By the definition of g(x), we have that 

g{yi + c/) 
= t? • (.9(y/ + c) -g{ci)) 



( 



\ 



i=o j=o j=o 



\ '/•-f.i 



■y' 



'4)1,2 



I 
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where Lp is either or a of polynomial whose least order is not smaller than 2. 
Thus, 

lim - — - = . 

lly/IKo ||y/|| 

Let X]iG7 ^i ' yi ~ '^7 ■ ("04,1 + ''Pi.2), and we obtain equation (|14l) in Proposition 
fT2l namely 



?(y/ + c/) « i7 • ^(V'i.i + V'i,2) 



ie/ 



Proof (Proof of Theorem\17\). It holds that for each i G / 
Kj = ■^(max(hi) - min(hi)) 
Ku, = -^ ^ w(i)(max(h,;) - min(hi)) 

where w{i) = /\i/Zi and A — X^jg/ '^j- Equation (PT)) follows. 



